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Abstract 

We present here a comprehensive search for the structure of the Si(553)-Au reconstruction. 
More than two hundred different trial structures have been studied using first-principles density- 
functional calculations with the SIESTA code. An iterative procedure, with a step-by-step increase 
of the accuracy and computational cost of the calculations, was used to allow for the study of this 
large number of configurations. We have considered reconstructions restricted to the topmost 
bilayer and studied two types: i) "flat" surface-bilayer models, where atoms at the topmost bilayer 
present different coordinations and registries with the underlying bulk, and ii) nine different models 
based on the substitution of a silicon atom by a gold atom in different positions of a vr-bonded 
chain reconstruction of the Si(553) surface. We have developed a compact notation that allows us 
to label and identify all these structures. This is very useful for the automatic generation of trial 
geometries and counting the number of inequivalent structures, i.e., structures having different 
bonding topologies. The most stable models are those that exhibit a honeycomb-chain structure at 
the step edge. One of our models (model f2) reproduces the main features of the room temperature 
photoemission and scanning-tunneling microscopy data. Thus we conclude that f2 structure is a 
good candidate for the high temperature structure of the Si(553)-Au surface. 

PACS numbers: 73.20.At, 71.15.Mb, 79.60.Jv, 68.35.Bs 
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I. INTRODUCTION 



During the last few decades several complex surface reconstructions that exhibit one- 
dimensional metal wires have been extensively studied. The interest in these systems has 
been mainly driven by the theoretical prediction that the behavior of a one-dimensiona 
metal should sensibly deviate from the Fermi liquid picture due to electron correlation. |l] 
However, the coupling between electrons and lattice vibrations is also enhanced in one di- 
mension and these systems present a strong tendency to suffer structural transitions that 
reduce their symmetry and destroy their metallicity. The competition between the effects 
of electron-electron and electron-phonon interactions makes the physics of one-dimensional 
systems quite a rich and interesting topic. In this context, the quasi-one-dimensional recon- 
structions formed by several metals on fiat and vicinal Si(lll) surfaces have been regarded 
as ideal systems to realize a one- dimensional metal. On the one hand, the existence of an 
electronic gap in the substrate allows for the existence of purely one-dimensional electronic 
states near the Fermi energy associated with the metal chains in the surface. On the other 
hand, it was believed that in some cases the rigidity of the substrate could allow for the 
observation of characteristic features of the Tomonaga-Luttinger liquid in the photoemis- 
sion experiments, [l, 2] avoiding the appearance of broken- symmetry phases at sufficiently 
low temperatures. Furthermore, ordered arrays of well aligned wires can be formed over 
large areas on these susbtrates, which allow for accurate measurements using angle resolved 
photoemission (ARPES). In the case of vicinal substrates the average terrace width can be 
controlled by the miscut angle, providing a route to tune inter-wire distances and interac- 
tions. ^ 

Some examples of quasi-one-dimensional reconstructions that have attracted much at- 



tention in recent years are the Si(lll)-(5x2VAu, 3 
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7| Si(557)-Au, fl, lo, Q, 
2^ and Si(lll)-(4x2)- 



371] surfaces. All these systems show 



clear one- dimensional features both in ARPES and scanning tunnelling microscopy (STM) 
experiments. At least in the cases of the Si(557)-Au, Si(553)-Au and Si(lll)-(4x2)-In sur- 
faces some of these one-dimensional bands seem to be metallic at room temperature. How- 
ever, these three reconstructions suffer structural distortions associated with metal-insulator 
transitions when the temperature is lowered. 
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The photoemission spectrum of the Si(557)-Au surface is dominated by two proximal one- 
dimensional bands with a parabolic-like dispersion. According to the ab initio calculations 
based on the existing structural model 10|, llll] these two bands are due to the spin-orbit 



splitting of a band coming from the hybridization of the states of the gold chains with their 
silicon neighbors. 13|, |38| These bands are half-filled and exhibit, as the temperature is 
lowered, a metal-insulator transition that is accompanied by a periodicity doubling along 
the steps in the STM images. [l2, Q, Q, Q 

The Si(553)-Au surface, with a smaller terrace width, shows two bands similar to those 
observed in the case of the Si(557)-Au system. It also shows a strongly-dispersing one- 
dimensional band with a fractional (~ |) filling, js], l3] It was proposed that this one-third- 
filled band could provide the appropriate conditions to observe large spin-charge separation. 
This is in contrast with half-filled bands, which are unstable against a Mott-Hubbard tran- 
sition for large values of the electron-electron interaction, preventing the observation of a 
Luttinger metal [3|]. Unfortunately, the Si(553)-Au surface also undergoes a metal-insulator 
transition at low temperatures. The Si(553)-Au surface exhibits two Peierls-like distortions 
with a x2 and x3 increase of the unit cell of the chains along the steps. This is consistent 
with the occupations of the bands mentioned above. Furthermore, the x3 periodicity seems 
to be associated with the step-edge, while the x2 periodicity appears in the middle of the 
terraces. This is consistent with the origin of the half- filled bands in the gold wires. Ab 
initio calculations indicate that the gold atoms prefer to occupy in silicon substitutional 
positions in the middle of the terraces. 

A good structural model is necessary to understand in detail the physics of this system. 
However, in contrast to the case of the Si(557)-Au surface, the structure of the Si(553)-Au 

n 

reconstruction is not well established. Recently Ghose et al. [19| proposed a structure based 
on the analysis of x-ray diffraction measurements. This model contains a double row of 
Au atoms decorating the step edge. This seems inconsistent with the experimental gold 

n n 

coverage and the model has proven to be unstable. [22i] In a recent work [18|] we studied, 
using first-principles density-functional calculations, five different structural models based 
on an analogy with the structure of the Si(557)-Au surface and an earlier proposal by Grain 
et al. [3l] with a single gold chain. The most stable structures were those in which the step 
edge is formed by honeycomb- chain. [3^ 

In the present paper we perform a more comprehensive search of the structure of the 
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Si(553)-Au reconstruction. We have considered models in which the topmost bilayer contains 
eight atoms per 5x1 unit cell, including one gold atom, with all different coordinations 
and registries with the underlying bulk. We refer to these models as "flat" surface-bilayer 
structures and we have found that there are 210 inequivalent structures of this type when 
reasonable physical constraints are applied. These structures are generated automatically 
and an iterative procedure, with a step-by-step increase of the accuracy and computational 
cost of the calculations, is used to study this large number of configurations. We have also 
considered nine different structures obtained by the substitution of a silicon atom by a gold 
atom in different positions of a 7r-bonded chain reconstruction of the Si(553) surface. 

As we will see, our most stable structures belong to the set of flat-bilayer models and ex- 
hibit a honeycomb-chain structure in the step edge. This in agreement with the observation 
of our previous study. Indeed the two most stable structures (f2 and f4 models) correspond 
to structures that have already been explored in Ref. Q. This confirms that the analogy 
with the Si(557)-Au surface provides a good route to study the structure of the Si(553)-Au 
reconstruction and, most probably, other related surfaces. [40|] The simulated STM images 



of the f2 and f4 models are in good agreement with the available experimental data. 



13, 



2l|41| 



m In 



Their band structures also reproduce the main features of the photoemission data, 
particular, the band structure of the f2 model shows two dispersive one-dimensional bands 
and one of them, associated with a silicon dangling-bond in the surface, has a small frac- 
tional occupation in agreement with the measurements. Therefore, we conclude that the f2 
model is a good candidate for the structure of the Si(553)-Au. 

Finally, we have developed a compact notation that allows us to label and identify all the 
structures. This notation is instrumental for the automatic generation of trial geometries 
and for counting the number of inequivalent (having different bonding topologies) structures. 
This notation is applicable to other similar surfaces and we think that it will be useful for 
other structural investigations. 
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II. METHODS 



A. Density functional calculations 



Most of the calculations were performed using the SIESTA code 
local density approximation 
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44j ] and norm conserving pseudopotentials j45|] . The gold pseu- 



43| We used the 



dopotential included scalar relativistic effects and was similar to that used in Ref. 



our previous calculations 



0,y 
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and in 



22|. Different basis sets were used for the silicon atoms 



see below for more details). For the fastest calculation a single-C basis (SZ) of numerical 



atomic orbitals 
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I48l | was used. The SZ basis uses only one radial function to represent 
the orbitals in the 3s shell and another one for the 3p shell. For more accurate calculations 
double-C basis (DZ), with two different radial functions for each angular momentum, and 
double-C polarized (DZP), including an additional polarization 3d shell, were used. For gold, 
a DZP basis was always used with two different radial functions to represent the 6s and 5d 
orbitals and a polarization 6p shell. The radius of the Si orbitals was 5.25 Bohr for those 
in the 3s shell and 6.43 Bohr for those in the 3p and 3d shells. For Au the radius of the 6s 
and 6p orbitals was 6.24 Bohr and 4.51 Bohr for the orbitals in the 5d shell. 

We modelled the surface using slabs with different thicknesses (depending on the desired 
accuracy). For the initial and fastest simulations we used a slab containing only two silicon 
bilayers, whereas four bilayers were used to more accurately explore the energetics of the 
most favorable structural models. The atoms in the bottom silicon layer were kept at their 
ideal bulk positions and saturated with hydrogen atoms. A vacuum gap of 15 A between 
neighboring slabs was used. To avoid artificial stresses the lateral lattice parameter was 
always fixed to the bulk theoretical value calculated with the similar approximations. The 
lattice parameter depends strongly on the basis set used. The calculated values are 5.53, 5.48, 
and 5.42 A for SZ, DZ and DZP basis sets, respectively. Therefore, it is necessary to rescale 
the coordinates when the result of a relaxation using a given basis set is used as the initial 
guess for a calculation with a more complete basis. Systems were relaxed until the maximum 
force component along any direction was less than 0.04 eV/A. Different Monkhorst-Pack 49| 
samplings of the surface Brillouin zone were utilized. They will be referred along this paper 
as MxN, where M refers the direction parallel and N perpendicular to the step edge. A 
real-space grid equivalent to a 100 Ry plane-wave cut-off was used. The simulated STM 
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images were obtained using the Tersoff-Haman theory 50| and the DZP basis set. 



51 



52|. 



For the most stable systems we also performed calculations using the VASP code 
Projected-augmented-wave potentials and a well converged plane-wave basis set with a cutoff 
of 312 eV were used. All structures were relaxed (the equilibrium lattice parameter of bulk 
silicon obtained with VASP was 5.41 A). 



B. Strategy of the structural search 

Searching for the lowest-energy configuration of any surface reconstruction is a formidable 
task. Almost any trial geometry that is not completely absurd relaxes to a nearby local min- 
imum. One can then generate manually new configurations from the most favorable ones 
using different physical (electron-counting, coordination, etc.) arguments, hoping to even- 
tually end up in the global minimum. However, the number of local minima scales roughly 
exponentially with the number of atoms involved in the reconstruction and the problem 
becomes intractable for large cells. Apart from this heuristic approach one can use more 
sophisticated algorithms that can automatically find optimum geometries. Some examples 



are simulated annealing 
genetic algorithms 
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|54| and Monte Carlo simulations of different types, including 
57|. Monte Carlo techniques have been traditionally used to 



find c^ 
tions 
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601 ]. Unfortunately these methods are computationally very expensive requiring 
long simulations with thousands of evaluations of the system energy (and interatomic forces 
in some cases). Particularly, genetic algorithms are very powerful but typically need hun- 



drends of generations, each one containing tens of trial geometries [58|]. Therefore, they are 
mostly restricted to the use of empirical interatomic potentials which, however, might not 
be sufficiently accurate to reproduce the energetics of the different geometries explored. 

In the present work we adopt a compromise between these two ideas. Due to the lack of 
reliable empirical potentials to represent the interaction between the gold and silicon atoms 
in the surface, we need to explore the energetics of the different models of the Si(553)-Au 
reconstruction at the density functional or similar level of theory. This precludes the use of 
Monte Carlo algorithms to perform a global search of the reconstruction structure. However, 
we do not want to restrict our search to explore a "few" structural models. Thus, we will 
rather make a systematic search of the optimum surface model within a large family of 
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physically motivated structures. Our approach is the following: 

i) a family of likely structural models for the Si(553)-Au surface reconstruction is defined 
using a heuristic procedure based on the analogy with other related and better known 
surfaces; 

a) a compact notation is designed to label unambigously each of the possible structures 
within this family; 

Hi) from each of these labels, a trial geometry is generated automatically and an initial 
constrained relaxation is performed to avoid the appearance of unphysical bond distances; 

iv) "fast" density functional calculations using SIESTA are used to relax each of the 
structures to its closest local energy minimum; 

v) the most stable configurations from step iv) are studied using more time consuming 
"accurate" SIESTA calculations; 

vi) finally, since the energy differences between different structural models are quite small, 
we check the energy ordering of our most stable configurations using a different methodology: 
plane-wave calculations using the VASP code. 



Dasis 



43|), 



By "fast" SIESTA calculations we mean here calculations performed using small 
sets (i.e, with a small number of basis orbitals per atom, like SZ or DZ basis sets 0, 
limiting the number of k-points and optimized degrees of freedom and/or using less stringent 
criteria than usual for the convergence of the self-consistency cycles. These approximations 
substatially reduce the computational cost. Therefore, it becomes possible to relax the 
hundreds of different configurations within our family of structures. We will see below 
that the energies obtained in this first step are reliable enough to select a set of a few 
tens of structures containing the most promising structural candidates. "Accurate" SIESTA 
calculations are performed for these configurations using DZP basis sets, a larger k-sampling 
and well converged self-consistency. We can see that the use of a code that utilizes basis 
sets of atomic orbitals is instrumental for this gradual increment of the accuracy of the 
calculations: while the pseudopotentials, density functional, and basic numerical scheme 
remains unchanged, the size of the basis set (the main factor limiting the size of the studied 
systems and the computational time) can be varied. This provides a very convenient way to 
deal with the trade-off between computational speed and accuracy of the calculation. For the 
Si(lll)-(5x2)-Au we have recently shown that SIESTA calculations using DZP basis sets are 
in excellent agreement with well converged plane-wave calculations, j?! In the present work 
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we confirm this observation for tlie proposed models for the Si(553)-Au surface. This seems 
to confirm that our final structures are the most favorable within the family considered here. 

In the following we present in detail the hierarchy of approximations used to perform 
our simulations. As a starting point we automatically generate approximate coordinates 
for all possible structures fulfilling certain conditions. These conditions will be explained in 
detail in the next section. Our family of structural models is based on plausible analogies 
with the structure of other related surface reconstructions like the Si(557)-Au or Si(lll)- 
(5x2)-Au. While reasonable atomic coordinates within the plane of the terrace (that we 
take as the xy-plane) are relatively easy to guess due to the registry with the sub-surface 
bilayer, the height of the different atoms in the surface bilayer is more problematic. For this 
reason, in the first relaxation (named Sz hereafter) only the atoms in the topmost bilayer 
are allowed to relax in the z-direction (normal to the terrace). This relaxation step ensures 
the interatomic distances to be reasonable without changing the topology of the surface 
bilayer and its registry with the underlying atoms. For the 5*2; relaxations we use a SZ basis 
set and a 2 x 1 k-sampling. To further accelarate the simulations the parameter determining 
the convergence of the density matrix in each relaxation step (DM. Tolerance 6l|]) is set to 
10^'^. The typical value given to this parameter to ensure a very good convergence of the self- 
consistent solution is 10~^. However, we have checked that increasing this value to 10~^ only 
introduces small errors in the calculation of energy and forces: for example, the maximum 
force difference during the relaxation of a few representative structures of those studied here 
was less than 0.01 eV/A when the two different convergence criteria were used. In spite of 
this moderate effect on the results, in some cases increasing the value of DM. Tolerance to 
10^'^ reduces considerably the number of iterations per self-consistency cycle. In the next 
step (named Sy hereafter) all atoms are allowed to move (except the silicion atoms in the 
bottom of the slab and the hydrogen atoms directly bonded to them). However, in order 
to preserve the topology of the selected configuration, the positions of the atoms along the 
direction parallel to the step edges (y-axis) is fixed. Other parameters have the same value as 
in the relaxation Sz. The purpose of the relaxations Sz and Sy is to provide a sound initial 
configuration from the coordinates generated automatically. Using this corrected guess we 
can proceed further allowing all degrees of freedom in the slab to relax and using a more 
complete DZ basis set and accurate 8x4 k-sampling. We call this the D* relaxation. Finally, 
in the D relaxations we further decrease the tolerance for the covergence of the elements of 
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the density matrix to its usual default value in SIESTA [61] of 10~^. We use a DZP basis 
set for our most accurate relaxations (DP hereafter). Adding a polarization shell with d 
symmetry can be especially important to accurately describe "unusual" coordinations of the 
silicon atoms which cannot be described with simple sp hybridizations. 

The use of this series of optimization schemes with ever-increasing accuracy 
{Sz^Sy-^D*^D^DP) is much more efficient than starting directly with a relaxation 
at the DP or similar level. The reason for this efficiency gain is twofold: i ) the initial 
coordinates used to start each relaxation have been optimized at the previous level and, 
therefore, they are an initial guess of increasing quality; ii) the energy estimates obtained 
with the less accurate relaxations, already at the Sy level, are accurate enough to allow 
discarding many of the possible configurations in favor of the most favorable models. This is 
the case even if at the Sy level we use a minimal basis set, a thin slab of two bilayers and we 
fix the atomic coordinate parallel to the step edge in order to preserve the selected bonding 
topology. One has to take into account that a minimal basis for Si contains only four or- 
bitals, a DZ basis eight orbitals, and a DZP basis thirteen orbitals. Thus the computational 
cost changes dramatically when changing the basis set size. Additional gains are obtained 
by using smaller k-samplings and reducing the number of step in each self-consistent cycle. 



C. Structural models: a labelling scheme 

The structural models that we considered here for the Si(553)-Au reconstruction are 

n 

based on the analogy with other similar reconstructions p| for which there is more structural 
information. In particular, the vicinal Si(557)-Au surface reconstruction has been extensively 



studied in recent years, and its structure has been established from X-ray studies |10[] and 
DFT calculations, j^, [ill, 13] Like in the case of the Si(553)-Au surface, the steps of the 
Si(557)-Au run along the [110] and each terrace contains a monatomic chain of gold atoms. 
The main difference is the orientation of the miscut angle, which is the opposite in both 
cases, js] The size of the terraces is also different, it is larger in the case of Si(557)-Au allowing 



for the presence of a row of silicon adatoms running parallel to the step edge |10l . Ill| which 
seems to be absent in the case of the Si(553)-Au surface. Therefore, it is clear that both 
structures can be quite different (particularly the step-edge configuration), and one has to 
be careful when trying to translate structural information from one system to the other. 
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Fortunately we also have a considerably amount of information about the reconstruction 
induced by submonolayer deposition of gold in flat Si(lll). The structure of the Si(lll)- 
(5x2)-Au surface, corresponding to a larger gold coverage than those commented above, was 
established experimentally by Marks and Plass [62] using high resolution electron microscopy 
and heavy-atom holography. Recent first-principles calculations have shown that the real 
structure is probably somewhat different from this original proposal . |63| and might 



even depend on the concentration of silicon adatoms present in the surface, p, 64] However, 
there are several common elements to most of the proposed models, which are also shared by 
the Si(557)-Au reconstruction. The most important features of the most stable structural 
models of these surfaces are: 

i) the reconstruction only involves the atoms in the topmost bilayer; 

a) the gold atoms occupy substitutional positions in the surface layer, which are much 
more favorable than adatom-like sites; 

in) positions of gold in the middle of terraces are favored over step-edge decoration; 

iv) frequent appearance of the so-called honeycomb chain (HC) 391] structure. 

Fig [1] (a) shows a schematic view of an unreconstructed Si (553) surface. Taking into 
account the points i) and ii) we only explore here reconstructions generated adding an 
additional bilayer ontop of this unreconstructed substrate. One of the silicon positions is 
replaced by a gold atom. Different registries with underlying bilayer are allowed, as well 
as, the presence of HC structures. Fig [T] (b), (c) and (d) show a few possible structures. 
Structure (b) recovers the unreconstructed surface. Structure (c) presents stacking fault in 
the middle of the terrace with the accompanying surface dislocations with under- and over- 
coordinated atoms (indicated by arrows in Fig. [1]). In panel (c) the surface bilayer contains 
a HC structure in the middle of the terrace. Notice that the HC reconstruction also creates 
a stacking fault towards the [112] direction (i.e., towards the inner part of the terrace). This 
stacking fault has to be corrected in order to connect with the bulk structure. Therefore, 
it is necessary to introduce a surface dislocation with over-coordinated atoms (marked with 
an arrow). 

The importance of the HC structure for the gold induced reconstructions on vicinal 



nasizea 



jy Crain et al. [3]. It is known to form the step-edge 



Si(lll) has been recently emj 
of the Si(557)-Au surface |9|, [lOl, 111 
proposals for the structure of the Si(lll)-(5x2)-Au reconstruction 5|, LZl]. The HC structure 



161 ]. It is also a key ingredient of the most recent 
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was initially proposed to explain the low coverage reconstructions induced by some alkali 
metals (Li, Na, K), Mg and Ag on Si(lll) [39|. We see that the HC structure involves two 
unit cells of the unreconstructed Si(lll) surface, with one atom removed from the top Si 
layer. This flattens the surface and removes surface stress. The two silicon atoms in the 
center of the HC structure have a double bond, which is further stabilized by its hybridization 
with the dangling-bond of atom immediately below. In the case of reconstructions induced 
by alkalis, the metal atoms donate electrons to the HC becoming a closed-shell structure 
and thus contributing to the stabilization of the surface. 39| In the case of gold, which has 
a stronger electron affinity, the situation is different. Gold is likely to take electrons away 
from the silicon structure. In principle, this does not prevent an electronic stabilization 
mechanism: one electron may be transferred to the 6s Au state, leading again to a closed- 
shell structure. However, we are typically far from this ionic situation. The states of gold are 
strongly hybridized with those of the neighboring silicon atoms creating several dispersive 
bands that are, in principle, metallic [sl, [ll, 1^ and dominate the photoemission spectra. 



In Ref. 



18l we explored a few structural models for the Si(553)-Au reconstruction based 



on an analogy with the Si(557)-Au surface. Here we want to move a step further and to 
make a comprehensive search among the structural models that can be built following the 
rules {i)-{iv) presented above. We consider all possible structures where the atoms of the 
topmost bilayer present coordinations between 2 and 4 with other atoms in the same bilayer. 
The final coordination depends on the registry with the underlying silicon structure. One of 
the silicon atoms in the unit cell is replaced by a gold atom. The Si(553)-Au reconstruction 
is known to suffer several distortions that increase the size of the unit cell along the step 



direction as the temperature is decreased. 2l|, |2J] However, here we only consider models 



that preserve the xl periodicity of the silicon substrate along the steps and, therefore, are 
relevant to model the high temperature structure. 

We have developed a simple labelling scheme for the family of structural models that 
fulfill the criteria presented above. We can label each structure and thus count the total 
number of different trial structures within this family. Furthermore, this scheme can be 
easily translated into a procedure to automatically generate the trial geometries. The basis 
of our labelling procedure can be found in Fig. [2l First, the possible positions within the 
xy plane of the surface atoms are discritized and approximated by the points of a grid. 
The grid is formed by nine columns and two rows. The nine columns correspond to the 
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positions of the atoms along the [112] direction in the terrace of an unreconstructed Si(553) 
surface (see the Fig. [1] (b)). Second, all possible structures created by distributing the 
atoms among the grid points can be translated into a sequence of nine numbers. The 
position along the horizontal coordinate (column) is indicated by the order in the numerical 
sequence. The first number corresponds to the atoms at the step edge. For a given column, 
a "2" ("4") indicates that a silicon (gold) atom is located in the higher row, i.e., in the 
middle of the rectangular terrace unit cell, whereas a"l" ("3") indicates that a silicon (gold) 
atom sits over a grid point in the lower row. A "0" indicates that there are no atoms in 
that column. Using this scheme, the unreconstructed Si(553) surface in Fig. [1] (b) can be 
label as (1,2,2,1,1,2,2,1,1), while structures in panels (c) and (d) would receive the label 
(1,2,2,2,1,1,2,1,1) and (1,2,2,1,0,1,2,1,1), respectively. Other examples, corresponding to 
models already studied in Ref. can be found in Fig. [2j 

In principle, using our notation we can generate M different models, with 

M = 2^^s.+i) ^ (1) 

{8-Ns.)\Ns.\ 

and Nsi being the number of silicon atoms in the surface bilayer. Since we always have one 
gold atom, the total number of atoms in the terrace unit cell is Natm=Ns'i+l. We consider 
here structures with Nsi=7- In this way, the family of structures studied here includes the five 

8 is a necessary condition 

to allow for the formation of the HC structure, which is one of the common building blocks 
to several gold induced reconstructions in Si(lll) and vicinal Si(lll) surfaces (see above). 
With Nsi=7 we have M=18432 different models. This large number can be considerably 
reduced imposing a few constraints to ensure that the models represent physically sound 
structures. These constraints are: a) in order to connect with the bulk structure the last 
number of the series must be either 1 or 3; b) the dangling-bonds of the underlying silicon 
bilayer must be saturated either by an atom or by a dimer as in the HC structure, therefore 
the first number of the labelling sequence must be always 1 or 3, c) the fifth number must 
be 0, 1 or 3 (if 0, then the neighboring numbers must be either 1 or 3), and d) the third and 
seventh number must be 0, 2 or 4 (if 0, then the neighboring numbers in the sequence must 
be 2 or 4); e) to ensure the connectivity within the surface bilayer, a non-zero number in the 
sequence cannot be surrounded by zeros and f) two or more zeros cannot appear together. 
Taken into account conditions (a)-(/), the number of possible configurations with Natm=^ 
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and N5j=7 is reduced to M=210. In this work we will explore these 210 configurations 
using "fast" SIESTA calculations, as decribed in the previous section, and a few tens of the 
most stable configurations will be selected to perform more accurate SIESTA and VASP 
calculations. 

Our notation provides information about the connectivity within the surface bilayer and 
the registry with the substrate. From a given sequence of nine numbers we can generate a 
trial geometry. However, we lack information about the heights of the different atoms. Due 
to this and to the discretization of positions in the xy-plane, the bond lenghts and angles in 
the automatically generated structures can considerably depart from the correct values. For 
this reason, as a first step to get sound initial configurations we need to perform constraint 
relaxations that, while preserving the bonding topology of the selected configuration avoid 
unphysical bond distances and angles. We use the Sz and Sy relaxations described in the 
previous section for this. 



Besides the family of structures described above, we have explored a 



ew structural models 



based on the 7r-bonded chain reconstruction of the Si(lll) surface. [65|, |66|, |67|] In principle, 
our notation cannot describe these bonding pattern: it can only describe structures which 
are based on a "fiat" surface bilayer. This is partially due to the lack of information about 
the atomic heights. However, we can modify our notation to describe the 7r-bonded chain 
structures. This is done allowing for a double occupation of the columns and is schematically 
illustrated in Fig. [3l These double occupation indicates the position of the vr-bonded chain 
in the structure. We still have the ambiguity about the relative height of atoms in the 
TT-bonded chain, which is known to be tilted. There are two possibilities which are usually 
referred as negative or positive tilt. 68| In our notation these two different tilts of the vr- 
chain are indicated by the order of the pair of indices, the second index corresponding to 
the higher atom. In Fig. [3] we present a Si(553) surface reconstructed with the negative 
tilt chain. This negative tilt yr-chain block corresponds to the label (...,1,0,2,21,1,...), while 
the label (...,1,0,2,12,1,...) denotes the positive tilt structure. Both configurations are quite 
similar and first-principles calculations predict them to be almost degenerate in energy and 



separated by a very small energybarrier. 69 



3 



Experimentally the positive tilt structure 



has been traditionally favored. [66|, l67|, l68| In the case of the Si(553) stepped surface, our 



calculations predict the negative tilt structure to be slightly more stable. 

This notation opens the possibility of generating and studying all possible structures 
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containing the vr-bonded chain. However, we have not pursued this approach here and we 
hmit to consider nine different structural models that are obtained after the substitution of 
one silicon atom by a gold atom in different positions of a vr-bonded chain reconstruction 
similar to that shown in Fig [31 

III. RESULTS 

We take advantage of the methodology described in sections IIIBI and III CI to make an 
extensive search of the structure of the Si(553)-Au reconstructions. We make a system- 
atic search within the 210 structures that can be generated with the notation presented in 
Sec. Ill CI (models based on "flat" bilayers with different coordinations and registries with 
the susbtrate) with seven silicon atoms and one gold atom in the terrace unit cell. We also 
present results from a much more restricted search for structures based on the vr-bonded 
chain reconstruction. Finally the most stable structures from these two searches are studied 
using accurate SIESTA and VASP calculations. We present results for the band structure 
and the simulated STM structures of some of the final models. 

A. Systematic search: "fiat" bilayers witii N„t,„=8 

We first explore the energy of the 210 possible configurations using our fastest relaxation 
schemes, Sz and Sy, described in Sec. Ill Bl These calculations transform an initial structure 
automatically generated from a given label into a physically sound structure. In spite of the 
thin slab and minimal basis utilized, the relative energies obtained at the Sy level (AEi) 
already provide a good guide to eliminate the most unstable structures. In Table [T] we can 
find a list with the 80 most stable configurations (AEi <33 meV/A^) obtained after Sy 
relaxations. Several of the initial structures converge to the same configuration, so these 80 
trial structures give rise to only 40 clearly different structures. This is seen in Fig.(|5]), where 
the plateaus in the energy curve correspond to this "lumping" of several initial geometries 
into a single geometry. This transformation typically takes place by a displacement of the 
surface bilayer as a whole, thus changing its registry with the underlying substrate, or by 
the movement of a vacancy to a neighboring position (using our notation this corresponds 
to a transformation (..,0,1,2,..) (..,1,0,2,..)). This happens for example in the case of 
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the (1,2,2,1,1,0,4,1,1) initial structure, that transforms into a configuration that is better 
described with the label (1,2,0,2,1,1,4,1,1) and is one of the most stable structures. This can 
be seen in Fig. H] (see also Tab. [I]). 

The 68 most stable structures, as predicted using AEi, are then calculated again, this 
time with more accurate relaxation schemes up to the DP level (for details, see Tab. [I]). 
Hereafter, we will refer to these final relative energies as AE2. to the DP level. The A 
comparison between AEi and AE2 can be found in Tab. [T] and Fig. [51 Of course, there are 
important differences between the relative energies calculated using the minimal basis and 
constrained relaxations {Sy scheme), and those obtained from the fully relaxed structures 
using DZP basis sets [DP scheme). However, the overall trend of increasing energy is 
recovered with the more accurate calculations. We also observe in this case that several 
initial geometries give rise to the same or very similar final structures. However, some of 
the approximate degeneracies found with the constrained relaxations are removed. 

Some of the most stable final geometries are marked with labels fl-flO both in Fig. [5] 
and Tab. [H These structures are shown in Fig. [6l The fast Sy relaxations (AEi) predicted 
fl to be most stable configuration. It exhibits a HC structure at the step-edge, while the 
gold atom is located at a surface dislocation in the middle of the terrace. The presence of 
a surface dislocation is necessary to recover the bulk stacking disrupted by the HC. The 
position of the gold atom seems reasonable, gold should be a better option than silicon to sit 
at the dislocation since gold does not exhibit strong directional bonding. However, using a 
more complete basis set (already a DZ basis gives the correct result) and a thicker slab the 
f2, f3 and f4 geometries become the most favorable structures (they are almost degenerate 
and ~4 meV/A^ more stable than fl). This points to the importance of using more complete 
(and thus fiexible) basis sets when the coordination of the surface atoms departs from simple 
sp hybridization (like in the case of the HC structure or at the surface dislocations). 

Our results indicate that the configurations featuring a HC structure at the step edge 
(similar to the Si(557)-Au structure) are the most stable, at least within the famil y o f 
reconstruction considered here. This confirms the results obtained in our previous work |l8| . 
where we only studied six different models (one of them was the model proposed by Crain et 
al. [3] that exhibits a x3 periodicity and, therefore, lies outside the present family of models). 



The configurations 
obtained in Ref. 



4 and f2 in the present study correspond with the most stable structures 
(named respectively I and II in that reference). Also configurations fS, 
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f9 and flO correspond to structures III, IV and V in Ref. isj. Configuration f3, however, is 
a new configuration exhibiting a double honeycomb structure similar to that found in some 
models of the Si(lll)-(5x2)-Au reconstruction. [5|, [3] 

In Sec. Ill CI we pointed out that our labelling scheme excludes, in principle, structures 
based on the 7r-bonded chain reconstruction. However, in Fig. we can find one structure 
where the 7r-bonded chain reconstruction has emerged spontaneously. In configuration f6 the 
gold atom is located very close to the step edge. The initial structure corresponds to a largely 



unreconstructed terrace. It is well known 



71| that the energy barrier for the transformation 



from the unreconstructed Si(lll) to the vr-bonded chain Si(lll)-(2x 1) reconstruction is very 
small. Therefore, the appearance of the 7r-bonded chain in this case is not very surprising. 

The experimental electronic band structure, as determined by photoemission experi- 
ments jsl, 21 1 , presents three bands with parabolic dispersion and strong one-dimensional 



character. Two of them are similar to those found for the Si(557)-Au surface and, therefore, 
can be assigned to the spin-split bands formed from the hybridization of Au 6p states with 
the sp lobes of the neighboring Si atoms, as proposed in our recent work. |l3] The third 
band appears centered around the same point in reciprocal space, but at higher energies and 
thus with a lower occupation around |. 

In Fig. [7] we can find the band structures of the models f2 and f3. The band structure of 
the model f4 has been published elsewhere [l8j] . The geometries of f2 and f4 are very similar 
(see Fig. [6]) but their band structures present some small but essential differences. Both 
models present one dispersive one- dimensional band coming from the hybridization of gold 



with its silicon neighbors. 



the Si(553)-Au surfaces [17, 



lis band can be identified with the spin-split bands observed in 
38 1 and Si(557)-Au js]. Two other surface bands appear close 



to the Fermi level: a dispersive band coming from a silicon dangling-bond in the surface and 
a band derived from the atoms at the step edge (that has a HC structure). These bands are 
also present for both models. However, while in the case of the f4 structure the step-edge 
band is completely filled and the dangling-bond band is empty [18:], in the case of the f2 
model these two bands cross and the dangling-bond band has a small occupation of ~ 
electrons closer to the experimental observation. 

The population of the different bands can be roughly understood taking into account the 
larger electron affinity of gold and the HC structure pS] as compared to other atoms in the 
surface. The population of the dangling-bond band is thus depleted in favor of the other 
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surface bands. The differences between the f2 and f4 structures are more subtle. These two 
models only differ in the position of the surface dislocation. Energetically this structural 
change has small consequences and both structures are almost degenerate. However, it 
changes the occupation of the dangling-bond band. In the case of the f2 model the dislocation 
involves the gold atoms and some of the silicon atoms of the HC structure. This increases 
slightly the energy of the bands associated with the HC and, as a consequence, the step-edge 
band transfers some of its population to the dangling-bond band. 

Fig. [7] also shows the band structure of the f3 model. The dangling-bond band is missing 
in this case. This spoils the comparison with experiment. The other two bands (gold and 
step-edge derived) are very similar to those found for the f2 and f4 models. This is reasonable 
taking into account the similar gold site and structure of the step-egde. 

Fig. [8] (b) shows the band structure of the f2 model calculated including the spin-orbit 
interaction with the VASP code. The band structure is in excellent agreement with that 
calculated using SIESTA. It also confirms that all the bands close to Ep with a significant 
weight in the gold atoms exhibit a splitting that has its origin in the spin-orbit interaction. 

Although the band structure of the model f2 does not exactly reproduce the photoemis- 

nnn 

sion results P, [171, |21[ , particularly the characteristic band fillings mentioned above, it has 
some clear qualitative similaritires with them. We find two bands with similar dispersions, 
with their minima ~1 eV and ~0.5 eV below the Fermi energy at the Brillouin-zone bound- 
ary. This corresponds very well to the photoemission data. Furthermore, the band with 
its minima at lower energy shows a notable band splitting near the Fermi energy due to 
the spin-orbit interaction. The is in good agreement with the experiments. However, the 
band structure of Fig. [8] (b) has two important differences compared to the photoemission 
data. Firstly, the dangling-bond band presents an important spin-orbit splitting associated 
with its appreciable hybridization with gold (see Fig. [7j). This splitting is not observed in 
the experiments. Secondly, the theoretical band structure has one extra band not seen in 
photoemission. This band is associated with the HC structure at the step-edge edge and 
has its minimum at F, contrary to the case of the other two bands. 

We can also compare the predictions for our models with the experimental STM images. 
In Fig. [9] and Fig. [10] we show the simulated STM images at different voltages for the f2 and 
f4 models. In general we can say that the comparison with the available experimental data 



is satisfactory. In agreement with the experimental images 2l|, |23|, l2J] the most prominent 
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feature is the step edge. Within the terrace we find signals coming from the row of gold atoms 
and its neighboring HC structure. The gold chain is seen as a continuous line for occupied 
states and presents more structure for empty states. For both positive and negative voltages 
we can also distinguish a signal coming from the unsaturated silicon dangling bond in the 
terrace close to the step edge. The step edge and the gold chain could be identified with the 
two parallel chains reported by Snijders et al. In their recent experiment, they observe 
a strong polarity despendence in the STM images, and in particular, find zigzag structures 
for empty states and ladder configurations for filled states. In figures [9] and [10] we have 
sketched some possible candidates for this kind of behaviour. For f2 (Fig. [9]), the ladder 
structure in filled states could result from the registry of the step-edge with respect to the 
gold atoms. For empty-states images the step-edge becomes slightly more visible and the 
zigzag geometry could result from the atoms inside the honeycomb-chain. In the case of 
model f4 (Fig. [TOj) . we can identify two entities that change their registries when going from 
filled to empty states. They are the step-edge and the complex formed by Au and the 
surface dislocation. They show, similar to the experiment, ladder and zigzag configurations 
respectively for filled and empty states. 



B. Restricted search: structures based on the vr-bonded chain 



We first explore two models of the Si (553) stepped silicon surface where the terraces are 
fully covered by a (2x1) yr-bonded chain reconstruction. They correspond to two slightly 
different arrangements of the vr-bonded chain (Fig [TT]) . Our most stable geometry (model 
pO) correspond to the so-called negatively tilted vr-bonded chain. {gsI] In low energy electron 
diffraction experiments of the flat _Si(lll) surface the positive-tilt 7r-bonded chain is usually 



favored over the negative tilt. 



66 



671 . |68| However, flrst-principles DFT calculations predict 



both structure to be ver y c lose in energy and there are conflicting claims about which of 



them is more stable. 



69 



7Q| 



We now proceed to make all possible substitutions of silicon by gold in the surface bilayer. 
This gives rise to nine different models for the Si(553)-Au reconstructions that we name pX, 
with X=l corresponding to a substitution at the step edge and X>1 to a substitution in 
the terrace. The flnal energies after accurate SIESTA relaxations {DP level) are listed in 
Tab. [TTl In several cases the initial structure was not stable and suffers strong modiflcations 
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during the relaxation. These changes are also summarized in Tab. |TT1 Model p4, illustrated 
in Fig. [12], is the most stable structure. The vr-bonded chain where the gold substitution 
takes place transforms into a structure similar to the unreconstructed Si(lll) surface. This 
is reasonable taking into account the atomic configuration of gold and confirms the tendency 
of gold to occupy substitutional positions in the middle of the terraces. 

m 



The next most favorable model, p2, is illustrated in Fig. [131 Again the vr-bonded chain 
where the substitution took place has disappeared. However, this time a configuration 
reminiscent of the HC structure has formed in the middle of the terrace. This transformation 
is accompanied by an expansion of the surface bilayer (see the forward movement of the step- 
edge atoms in Fig. [T^ . 

The band structure calculated for the p2 model, shown in Fig. [12] (c), presents similar 
characteristics to the experimental band structure (see the previous section). We can see in 
Fig. [13] (c) that the p2 model presents two metallic dispersive bands centered at point K in 
the zone boundary. One of theses bands is associated with the Si-Au bonds between gold 
and the neighboring silicon atoms in the step edge. It has a considerable gold weight and, 
therefore, is expected to exhibit a splitting if the spin-orbit coupling is taken into account. 
The other band, however, is mainly derived from the unsaturated dangling bonds of the 
silicon atoms at the step edge and presents a smaller filling. As expected, the vr-bonded 
chain structure that remains in the terrace does not give rise to any metallic band. 

The electron pocket of the step-edge band around K has an occupation of ~l/4, quite 
close to that found in the experiment. There is another small electron pocket associated 
with the step-edge around F which is not observed in the experiment. Thus we can assign 
a population of ~0.4 to the surface bands with a larger weight in the step-edge atoms. The 
dispersive band with a mixed silicon-gold character has an occupation of ~0.6 electrons. 
Therefore, the total population associated with the step-edge derived surface bands is one. 
This is somewhat surprising if we take into account that the Si atoms at the step-edge have, 
in principle, three electrons to populate these surface bands. However, one of these electrons 
is transferred to the Au 6s states (which appear several eV below Ep) and the other electron 
populates states with large contribution from the HC structure (which does not exhibit any 
metallic band). 

Further investigation of the p2 model reveals that it is a metastable configuration. Per- 
forming the relaxations with a more stringent force tolerance of 0.01 eV/A results in a more 
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stable structure. This structure, labelled p2* and illustrated in Fig. [HI exhibits a strong 
rebonding of the step edge. The vr-bonded chain suffers a small translation along the [110] 
direction in order to saturate the dangling bonds at the step edge. The new position of the 
atoms of the vr-bonded chain seems to be a compromise between creating a surface dislo- 
cation and saturating the dangling bonds at the step edge. The band structure, shown in 
Fig. [m (c), is quite similar to that of the p2 model. The role of the step-edge atoms is 
now played by the silicon atoms that formed the vr-bonded chain. A quite dispersive surface 
band, coming from these atoms, appears centered at K. Another dispersive band with a 
strong gold weight appears at lower energies. Although the topology of the band structure 
resembles that observed in the experiment, the filling of the parabolic silicon band is in this 
case close 0.4, i.e. larger than that observed experimentally. Fig. [T3] (c) also highlights one 
of the bands associated with the HC structure, showing a characteristic dispersion. 



The band structure of model p2*, calculated using VASP and including the spin-orbit 
interaction (Fig. [H] (a) ) shows quite good, although not perfect agreement, with the ex- 
periments. As in the case of the f2 model (see above), due to the spin-orbit interaction 
the dispersive bands suffer a splitting proportional to their weight in the gold atoms in the 
surface. Thus, the splitting is much larger for the dispersive band starting at lower energies. 
Again, in contradiction with the experimental observations we find some degree of splitting 
also for the parabolic band appearing at higher energies, although this splitting is smaller. 
The overall conclusion from Fig. [8] (a) is that the p2* outperforms the f2 model in terms of 
reproducing the photoemission results. 

The simulated STM images for this structural model are shown in Fig. [151 The depen- 
dence on the polarity seems to be much larger than for models f2 and f4. The structure 
of the STM images becomes more complex when going from filled to empty states as sev- 
eral "spots" appear and the identification of zigzag or ladder-like structures becomes quite 
arbitrary. 

Models p4 and p5 exhibit a similar rebonding of the step edge as explained above in the 
case of p2 and p2*. However, these new geometries, p4* and p5*, are not as stable as p2* 
(see Tab. [IT])- Geometry p5* also develops the HC structure in the terrace as can be seen in 




Fig.[Ii 
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C. Most stable structures: combined SIESTA and VASP results 



In Tab. IIIII we compare the converged energies of the most stable models found in the 
previous sections for the Si(553)-Au reconstruction. As discussed in the previous sections, 
these models have been found i) using a systematic search among all possible model based 
on a flat surface bilayer and eight atoms in the terrace unit-cell and, ii) the substitution of 
gold in different positions of a yr-bonded reconstruction of the Si(553) terraces. These two 
classes of models have different number of atoms. In order to compare the relative surface 
energies we need to define the chemical potential of silicon. Since the surface should be 
in equilibrium with bulk, we have chosen the chemical potential equal to the total energy 
of a silicon atom in bulk. Since the energy differences are quite small we have decided to 
perform the calculations of the most stable structure with a different methodology in order 



to cross-check our results. We have used the plane- wave code VASP for this purpose. |5ll.l52| 
We can see that there is an excellent agreement between SIESTA and VASP results. This 
was also observed in otKe. ..Ua. .westigation. Geo.et*. f4 and f2 ,ee Fi, O a.e 
the most stable structural models of those found in this extensive structural search. These 
two models can be considered degenerate within the precision of the calculations. 

Models f2 and f4 were already obtained as the most stable ones in a much more restricted 



structural search. 



18j The present calculations confirm that they are certainly among the 



most stable reconstructions of Si(553)-Au surface that only involve the topmost bilayer. 



IV. CONCLUSIONS 



We have presented a comprehensive study of the structure of the Si(553)-Au reconstruc- 
tion. We have considered reconstructions restricted to the topmost bilayer and studied two 
possible types: (a) "flat" surface-bilayer models, where the atoms at the topmost bilayer 
present different coordinations and different registries with the underlying bulk, and (b) nine 
different models based on substitutions silicon by gold in different positions of a vr-bonded 
chain reconstruction of the Si(553) surface. 

Although the unreconstructed Si(553) surface has nine inequivalent atomic positions, 
within the familiy of the "flat" surface-bilayer models, we have focused on structures con- 
taining only eight atoms (seven silicon and one gold atoms) in the topmost bilayer. We found 
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in previous investigations that this was crucial to reduce the number of danghng bonds in 
the final structures and to facilitate the appearance of the so-called h oney comb-chain struc- 
ture 9|]. The honeycomb-chain has been proposed by us I 111 . Ilq . llSl ] and by other 
authors 3|, |5| as one of the key ingredients of the reconstructions induced by gold on Si(lll) 
vicinal surfaces and the results presented here seem to confirm this idea. 

We have developed a compact notation that allows us to label and identify all the struc- 
tures. This notation is instrumental for the automatic generation of trial geometries and 
for counting the number of inequivalent structures, i.e., having different bonding topologies. 
There are several thousands of fiat surface-bilayers models with eight inequivalent positions 
in the topmost bilayer, however, applying a few physically motivated constraints we obtained 
only 210 different models. All these structures, along with those based on a 7r-bonded re- 
construction of the terraces of the Si(553) surface, have been studied using first-principles 
density functional calculations with the SIESTA code. An iterative procedure, with a step- 
by-step increase of the accuracy and computational cost of the calculations, was used to 
allow for the study of this large number of configurations. 

Some conclusions can be drawn from our investigation: 

i) Among all the explored models, the most stable configurations are those that present 
a honeycomb-chain structure at the step edge. In particular, models f2 and f4 are the most 
stable (see Fig. [6]) with almost degenerate total energies. The structure of both models is 
very similar, the sole difference steming from the relative position of the row of gold atoms 
and the surface dislocation induced by the presence of the honeycomb- chain. 

ii) Within this group of models, the band structure of model f2 shows a reasonable 
agreement with the photoemission measurements. 17|, |2l| . It exhibits two one-dimensional 
bands that share some similarities with the experimental results: both bands have a parabolic 
dispersion with the band minima at the K-point, one of them having a small fractional 
occupation (~ |). However, both bands present some gold character and thus exhibit an 
appreciable spin-orbit splitting. In the experiment only one band shows a clear splitting. 

in) The STM images and their dependence on the bias voltage is satisfactorily repro- 
duced by the simulated STM images of models f2 and f4. This is also in agreement with 
the observations of Ryang and coworkers |41i]. These authors have recently found a good 
agreement between the experimental STM images obtained for the Si(553)-Au surface at 
low defect concentration and those simulated for a model similar to our f4 structure. 



23 



iv) Of those models based on a 7r-bonded reconstruction of the Si(553) surface, only p2 
and p2* (p2* being energetically more favorable) present a band structure in reasonable 
agreement with the photoemission measurements. In particular, they exhibit two dispersive 
one-dimensional bands similar to the experimental ones: one of them has a filling close to ^ 
with a strong gold character, and thus exhibits a strong spin-orbit splitting; the other one, 
coming from the silicon step edge, has a smaller occupation and a quite small spin-orbit 
splitting in good agreement with the photoemission data. 

v) Simulated STM images of the p2* model present a strong polarity dependence. Un- 
fortunately, patterns appearing are hard to identify with those seen in the experiments. 

In summary, it seems that model f2 is a good candidate for the structure of the Si(553)-Au 
surface. This model presents the lowest total energy and its band structure and simulated 
STM images are in general agreement with the existing experimental data. The agreement, 
however, is not perfect. For example, model p2* seems to provide a better comparison with 
the photoemission data. Prom a more general perspective, our data clearly point to the 
importance of the honeycomb-chain structure to stabilize the surface and, probably, to the 
presence of a complex rebonding of the step edge, like that found in our pX* models. If this 
the case, the structure of the step edge could appreaciably depart from that of the initial 
models considered in the present study. More theoretical and experimental work is therefore 
necessary to determine if model f2 indeed corresponds to the high temperature structure of 
the Si(553)-Au surface reconstruction. 
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No. 


Initial 


AEi 


AE2 






configuration 


(meV/A2) 


(meV/A2) 


1 


1,2,2,0,1,1,4,1,1 


0.00 


0.00 




2 


1,2,2,1,0,1,4,1,1 


0.00 


0.00 




3 


1,2,2,1,1,0,4,1,1 


0.00 


0.00 




4 


1,2,0,2,1,1,4,1,1 


0.00 


0.00 




5 


1,0,2,2,1,1,4,1,1 


0.01 


-0.01 




6 


1,2,2,0,1,3,2,1,1 


5.51 


-4.27 




7 


1,2,2,1,0,3,2,1,1 


5.51 


-4.29 




8 


1,0,2,2,1,3,2,1,1 


5.51 


-4.26 




9 


1,2,2,1,3,0,2,1,1 


5.53 


-4.30 




10 


1,2,0,2,1,3,2,1,1 


5.54 


-4.28 




11 


1,2,2,1,1,0,2,3,1 


6.08 


-4.02 




12 


1,2,0,2,1,1,2,3,1 


6.09 


-4.01 




13 


1,2,2,0,1,1,2,3,1 


6.09 


-4.02 




14 


1,2,2,1,0,1,2,3,1 


6.09 


-4.01 




15 


1,0,2,2,1,1,2,3,1 


6.12 


-4.02 




16 


1,0,2,2,1,4,2,1,1 


6.21 


-4.42 




17 


1,2,0,2,1,4,2,1,1 


6.23 


-4.46 




18 


1,2,2,0,1,4,2,1,1 


6.29 


-4.45 




19 


1,2,0,2,3,2,2,1,1 


7.47 


1.75 




20 


1,2,2,0,3,2,2,1,1 


7.49 


1.73 




21 


1,0,2,2,3,2,2,1,1 


7.49 


1.74 




22 


1,0,4,1,1,2,2,1,1 


14.69 


9.55 




23 


1,2,2,1,1,0,2,4,1 


14.76 


-0.20 




24 


1,2,2,1,0,1,2,4,1 


14.76 


-0.20 




25 


1,2,2,1,1,2,0,4,1 


14.77 


-0.21 




26 


1,2,0,2,1,1,2,4,1 


14.77 


-0.17 


(£7) 


27 


1,2,2,0,1,1,2,4,1 


14.78 


-0.20 




98 


199119/ini 
i,Z,Z,i,i,Z,4l,U,l 


1 /I 


-0.18 




29 


1,0,2,2,1,1,2,4,1 


15.41 


-0.20 




30 


1,0,2,2,1,3,2,2,1 


17.49 


16.12 


28 


31 


1,2,0,2,1,3,2,2,1 


17.49 


16.13 





No. 


Initial 


AEi 


AE2 






configuration 


(meV/A2) 


(meV/A2) 


41 


1,2,2,0,3,1,2,1, 


E3.88 


7.47 




42 


1,2,2,3,0,1,2,1, 


E3.90 


7.48 




43 


1,2,2,3,1,0,2,1, 


E3.91 


7.47 




44 


1,2,0,2,3,1,2,1, 


E3.91 


7.47 




45 


1,0,2,2,3,1,2,1, 


E3.91 


7.47 




46 


1,2,0,2,1,1,2,2, 


24.13 


12.18 




47 


1,2,0,2,1,2,4,1, 


E5.23 


7.76 




48 


1,0,2,2,1,2,4,1, 


E5.24 


7.77 




49 


1,2,2,0,1,2,4,1, 


E5.24 


7.77 




50 


1,2,0,2,1,2,2,3, 


E5.41 


4.17 


(f8) 


51 


1,2,2,0,1,2,2,3, 


E5.41 


4.13 




52 


1,0,2,2,1,2,2,3, 


E5.42 


4.13 




53 


1,2,0,4,1,2,2,1, 


E5.67 


11.56 




54 


1,2,4,0,1,2,2,1, 


E5.72 


11.50 




55 


1,4,2,0,1,1,2,1, 


E6.31 


9.40 




56 


1,4,2,1,0,1,2,1, 


E6.34 


9.40 




57 


1,4,2,1,1,0,2,1, 


E6.36 


9.40 


(i9) 


58 


1,4,0,2,1,2,2,1, 


E7.35 


12.50 




59 


1,0,4,2,1,2,2,1, 


E7.35 


12.49 




60 


1,2,4,0,1,1,2,1, 


E7.37 


15.34 




61 


1,4,2,0,1,2,2,1, 


E7.40 


12.50 




62 


1,2,2,3,1,2,2,0, 


E7.62 


2.66 




63 


1,2,2,3,1,2,0,2, 


E7.63 


2.66 


(flO) 


64 


1,2,2,3,0,1,2,2, 


E7.63 


2.68 




65 


1,4,2,1,1,0,2,2, 


E8.08 


11.52 




66 


1,0,4,2,1,1,2,2, 


E8.09 


9.70 




67 


1,4,2,1,0,1,2,2, 


E8.09 


9.75 




DO 


1/1901199 


198 no 
jzo.uy 


9.75 




69 


1,2,2,1,0,3,2,2, 


E8.91 






70 


1,2,0,4,1,1,2,1, 


E9.48 






71 


1,2,4,1,0,1,2,1, 


E9.48 







Name 


Configuration 


AE (meV/A2) 


P 


(1,0,2,12,1,0,2,12,1; 




1.82 


pO 


(1,0,2,21,1,0,2,21,1; 




0.00 


Pl 


(3 2 21 1 2 21 r 




14.92 


p2 


(1,0,4,21,1,0,2,21,1' 


^(1,4,2,1,1,0,2,21,1) 


4.13 


p3 


(1,0,2,41,1,0,2,21,1' 


^(1,2,4,1,1,0,2,21,1) 


8.15 


p4 


fl 2 23 1 2 21 r 


^fi 2 2 3 1 2 21 1) 


0.00 


p5 


(1,0,2,21,3,0,2,21, r 


^(1,2,0,2,1,3,2,21,1) 


7.07 


p6 


(1,0,2,21,1,0,4,21,1; 


^(1,2,2,1,1,4,2,1,1) 


9.13 


p7 


(1,0,2,21,1,0,2,41,1; 




8.30 


p8 


(1,0,2,21,1,0,2,23,1; 




11.98 


p9 


(1,0,2,21,1,0,2,21,3; 




15.90 


p2^ 






-4.28 


p4'^ 






-2.69 


p5* 






0.38 



TABLE II: Relative surface energies of different structures based on a 7r-bonded chain reconstruc- 
tion of the terraces of the Si(553) and Si(553)-Au surfaces. Models p and pO correspond to the 
clean silicon surface. pX corresponds to the substitution of a gold atom in position "X" of structure 
pO. During the relaxation process several of these structures transform into configurations with a 
different bonding topology. This change is also indicated. pX* refers to pX configurations after 
the formation of bonds between the step edge and the neighboring 7r-chain structure. All energies 
correspond to our most accurate SIESTA calculations. 
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Name 


AE (meV/A2) 




SIESTA 


VASP 


p2* 


4.85 


4.93 


p4* 


6.44 


6.54 


p5* 


9.51 


9.13 


fl 


4.63 


4.27 


f2 


0.17 


-0.08 


f3 


0.42 


0.51 


f4 


0.00 


0.00 



TABLE III: Relative surface energies of our most stable models calculated using both SIESTA 
and VASP. 
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FIG. 1: (a) Schematic view of the unreconstructed Si(553) surface. The steps run parallel to 
the [110] direction and are oriented towards the [112] direction. The rectangle indicates the "unit 
cell" within the terrace, containing nine inequivalente silicon sites (with four unsaturated dangling- 
bonds). Panels (b), (c) and (d) show configurations generated by adding a silicon bilayer on top 
of the structure in panel (a). The different reconstructions explored in this paper are obtained by 
changing the structure of this surface bilayer and/or its registry with the underlying atoms. Open 
triangles represent the topmost atoms and solid circles the higher atoms in the second bilayer. 
These structures (and all the structures considered in the present work) preserve a xl periodicity 
along the step. Structure (b) recovers an unreconstructed silicon structure. Structure (c) presents 
a surface dislocation close to the step egde that generates a stacking fault (SF) that is later 
corrected creating another surface dislocation in order to connect with the bulk structure. Surface 
dislocations create under- and over-coordinated atoms which appear indicated by arrows. Model 
(d) presents a honeycomb-chain (HC) structure and the accompanying surface dislocation. 



FIG. 2: (a) and (c) show two possible structures for the Si(553)-Au surface already explored in 



Ref. Il8l . These structures are characterized by the presence of a HC structure close to the step- 
edge. Solid lines indicate the bonds between atoms in the topmost bilayer of a given terrace, thin 
dashed lines correspond to the underlying silicon bilayer, and thick dashed lines indicate a few 
bonds of the upper terrace. The small solid circles mark the positions of the higher silicon atoms 
in the underlying bilayer. The large open circles mark the substitutional sites occupied by the gold 
atoms. The rectangles mark the terrace "unit cell" . Panels (b) and (d) schematically explain how 
the structure of the surface bilayer can be translated into a sequence of nine numbers. First the 
possible positions of the atoms are approximated by the points of a grid. The grid is formed by 
nine columns and two rows. The position along the horizontal coordinate is indicated by the order 
in the numerical sequence. The first number corresponds to the atoms at the step edge. A "2" 
("4") indicates that a silicon (gold) atom is located in the higher row, i.e., in the middle of the 
rectangular cell. A "1" ("3") indicates that a silicon (gold) atom sits over a grid point in the lower 
row. A "0" indicates that there are no atoms in that column. 
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FIG. 3: (a) Scheme of a vr-bonded chain reconstruction of the Si(553) surface, and (b) the proposed 
notation for such structure. The presence of the vr-bonded chain at a certain location is indicated 
by the double occupation of the corresponding column, the second atom of this pair occupies the 
higher position along the vr-chain. 

FIG. 4: Relaxation of the structure generated from the label (1,2,2,1,1,0,4,1,1) at the Sy level. 
Insets (a) and (b) present a top and a lateral view of the automatically generated structure. Inset 
(c) presents the structure after the relative heights of the atoms have been corrected by the Sz 
relaxation. In point (d) we have already reached a structure similar to the configuration that we 
want to explore. However, this configuration is not stable and transforms by the displacement of 
the surface bilayer (inset (e)) into a different structure. The final structure, shown in inset (f), is 
better described by the (1,2,0,2,1,1,4,1,1) label. 

FIG. 5: Relative energies after fast Sy (AEi, solid circles) and accurate DP (AE2, open diamonds) 
relaxations for the systems listed in Tab. [H The system labels fl-flO correspond to those of Tab. [J 
and Fig. 

FIG. 6: Final geometries (at the DP level) for a few selected configurations from those listed in 
Tab. H] and Fig. [5j In the case of the f2 and f3 structures some symbols are assigned to a few atoms 
in the surface. They will be used to indicate the main character of the different bands (see Fig. [7|). 

FIG. 7: Band structures of models f2 and f3 calculated using SIESTA. The atomic character of 
different bands is indicated using the symbols in Fig. [6l Filled squares indicate the contribution 
coming from the gold atoms and their silicon neighbors, open circles that coming from the atoms at 
the step edge, and open diamonds that of some silicon atoms in the surface presenting unsaturated 
bonds. 

FIG. 8: Band structures calculated using VASP and including the spin-orbit coupling for (a) 
model p2* and (b) model f2. 
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FIG. 9: Simulated STM images, calculated using SIESTA, for model f2. The bias voltage is 
indicated in Volts in each panel. Some features having a strong polarity dependence are indicated 
by filled circles. 

FIG. 10: Simulated STM images, calculated using SIESTA, for model f4. The bias voltage is 
indicated in Volts in each panel. Some features having a strong polarity dependence are indicated 
by filled circles. 

FIG. 11: Two possible models of the Si(553) surface based on a 7r-bonded chain reconstruction of 
the terraces. Model p in panel (a) and model pO in panel (b). 

FIG. 12: Top (a) and side (b) view of the relaxed structure of the p4 model. 

FIG. 13: Top (a) and side (b) view of the relaxed structure of the p2 model, (c) Electronic band 
structure. The symbols highlight those surface bands with an appreciable weight from the surface 
atoms marked with the same symbols in panel (b). The inset shows a schematic representation of 
the most prominent surface bands. The HC structure in the middle of the terrace is indicated by 
a box. 

FIG. 14: Top (a) and side (b) view of the relaxed structure of the p2* model, obtained after the 
reconstruction of the step edge, (c) Electronic band structure. The symbols indicate those surface 
bands with an appreciable weight from the surface atoms marked with the same symbols in panel 
(b). The inset shows a schematic representation of the most prominent surface bands. The HC 
structure is highlighted by a box. 

FIG. 15: Simulated STM images, calculated using SIESTA, for model p2*. The bias voltage is 
indicated in Volts in each panel. 

FIG. 16: Side view of the p5* model. The HC structure emerges and is highlighted by a box. 
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